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Abstract A numerical model that predicts the stochastic 
dispersion associated with industrial kiln drying of timber 
was adapted to conventional drying and evaluated with ex¬ 
perimental data. The theoretical aspects of the model are 
briefly explained, a selection of the calibration parameters 
was carried out, and a new empirical dispersion factor is 
proposed to account for all unknown sources of random 
behavior. The model was calibrated with six experimental 
runs of western hemlock and amabilis flr (116 mm 2 timbers) 
to an average moisture content (target) of 14%-20%. It was 
found that with implementation of the dispersion factor, the 
number of required simulations is considerably reduced, the 
calibration results are consistent for all the experimental 
runs, and the target moisture content along with its standard 
deviation can be well reproduced using the all-run average 
parameters. 

Key words Conventional wood drying ■ Numeric stochastic 
model • Moisture dispersion prediction 


Introduction 

Conventional drying of wood is a process of random charac¬ 
teristics. In one commercial run, thousands of pieces are 
dried simultaneously, all with a varied degree of difference 
in their thermophysical properties and initial moisture 
content (/V/,). Furthermore, ambient air conditions are not 
completely homogeneous throughout the kiln where tem¬ 
perature and relative humidity usually oscillate in time and 
space due to control system limitations. As a result, timber 
positioning inside the kiln also affects the drying rate, and 
the surface or internal checks that appear in some or all of 
the timbers could change their internal heat and mass trans¬ 
fer characteristics. 
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As a result of these random factors, the final moisture 
content (M f ) is not an exact value but a distribution within 
a certain range of dispersion. The understanding and pre¬ 
diction of this intrinsic characteristic of the drying process is 
a real challenge to both structural timber producers and 
users, as the M f range of dispersion is maintained within 
specific limits (13%-19%); otherwise, thick timbers could 
exhibit unacceptable degradation (related to over-dried 
timber; i.e., M t < 13%) or structural instability (related to 
wet timber; i.e., M f > 19%). Different sort-out procedures 
have been proposed and tested to reduce the M t range of 
dispersion, 1 ' 2 but they require appropriate optimization to 
ensure that revenues due to the improved quality of the 
timber will cover the costs on the extra labor and equipment 
requirements. Therefore, an accurate stochastic prediction 
model would be useful for the optimization of these timber¬ 
drying strategies and the development of new ones. 

Although some other methods have been proposed to 
simulate stochastic drying of timber, such as the Monte 
Carlo model proposed by Kayihan, 3 they are based on a 
piece-by-piece random simulation, so they produce a 
slightly different result each time the model is applied to the 
same initial conditions. Because of this, the results with 
these models are also uncertain, so there is no easy way they 
can be implemented in an optimization procedure. There¬ 
fore, the stochastic model described in Elustondo et al. 4,5 
was adapted to simulate conventional drying of timber. The 
model is briefly described and then considerably simplified 
by introducing of a new empirical dispersion parameter. 
Finally, it is calibrated with experimental data and used to 
predict conventional drying of timber for a process in which 
the targeted M f is changed arbitrarily. 


Theory 

Stochastic model 

It is generally accepted in the wood-drying research com¬ 
munity that the Monte Carlo method 3 is a standard proce- 
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dure to analyze the outcome of stochastic processes. 6 The 
goal of this method is to generate a piece-by-piece random 
simulation by considering different locations inside the kiln 
and different timber thermophysical properties. Stochastic 
behavior is generated by a computer based on a random- 
number generator, and Weibull and Gauss probability 
distribution curves are used to describe M i and model pa¬ 
rameter dispersions, respectively. Although realistic, this 
approach produces a different result each time a new run is 
performed. The method is equivalent to the measurement 
of a sample containing a limited number of timbers, 250 in 
the case of Kayihan. 3 For this study, the stochastic model is 
developed to simulate not a particular experiment but the 
most probable result as if the data were averaged for an 
infinite number of runs. The proposed method is not far 
different from the Monte Carlo one; but instead of repeat¬ 
ing the simulations several times using random sequences of 
parameter values, the method takes samples from all the 
probability space to represent totally the possible param¬ 
eter combinations. 

The theoretical foundation of the method is explained by 
Elustondo et al. 2 Basically, the multidimensional space de¬ 
scribed by the model parameters is divided into a set of 
subspaces associated with a parameter combination (nor¬ 
mally the center point of the subspace) and a probability to 
occur. All generated parameter combinations are concur¬ 
rently simulated, and the resulting M f versus probability 
points are combined to obtain the M f distribution. Green M i 
distribution, another stochastic parameter, is treated differ¬ 
ently by the method to increase the efficiency of the simula¬ 
tion. As explained by Elustondo and Avramidis, 5 the M i 
distribution is also divided into a set of M i subintervals not 
associated with the center point but to the limits of each 
subinterval. 

The number of simulations required for this method is 
equal to the product of the number of points used to de¬ 
scribe each stochastic parameter. In practice, acceptable 
histograms can be obtained using five points for the timber 
properties and six points for the green M i distribution. This 
means that for a two-parameter model (such as the one 
proposed by Kayihan 3 ) 150 simulations are required. Al¬ 
though more points could be used in the representation of 
the stochastic parameters, the model considers all possible 
combinations, so the average and range of dispersion of the 
results is always more or less well represented. 


Intrinsic dispersion 

Conventional drying is affected by the differences in the 
timber characteristics, the kiln nonhomogeneous drying 
environment, the instability in the kiln-controlled variables, 
the differences in the sticker channels among successive 
timber packs, and the timbers with three faces instead of 
two exposed to the air, to name a few. For this reason, the 
stochastic model was complemented with an additional in¬ 
trinsic dispersion parameter representing all unknown 
sources of random behavior. This unknown source of dis¬ 
persion is applied to each simulated M f before the M f histo¬ 


gram is generated, so it does not increase the number 
of parallel simulations. For this particular study, the stan¬ 
dard deviation of the intrinsic dispersion ( d ) is assumed 
to be directly proportional to the amount of evaporated 
water. 

a = S a {M, - M) (1) 

Even though this is an empirical relation (which proved to 
be enough to predict M data dispersion under the condi¬ 
tions of this study), Eq. 1 can be derived by assuming that 
the drying rate range of dispersion on limbers starting with 
a similar M i is proportional to their average drying rate. 


Deterministic conventional model 


The described method incorporates stochastic behavior 
through the parameters of a deterministic model, but selec¬ 
tion of the deterministic model, simple or complicated, is a 
question of convenience. In the case of this study, the tim¬ 
ber model is based on well-known heat and mass transfer 
differential balances for porous materials. 7 9 The model is a 
set of three differential equations taking into account heat, 
vapor, and air balances; and for kiln drying it is solved with 
regard to the internal moisture content (M), temperature, 
and pressure profiles. Heat flux is replaced by Newton’s law, 
bound water and binary vapor-air diffusions are replaced by 
Fick’s law, capillary water and total gas fluxes are replaced 
by Darcy’s law, vapor equilibrium and capillary gradient 
pressures are replaced by their partial derivatives with re¬ 
spect to the temperature and M, and the ideal gas law is 
assumed for the gas phase. Finally, and after mathematical 
rearrangement and some minor simplifications, the follow¬ 
ing set of differential equations is obtained. 

^ = V-VM+ V-A^VEF V-A XP VF g (2) 

at 

~\nr 
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where the values of the coefficients are, above the fiber 
saturation point: 
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and below the fiber saturation point: 
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that conventional drying is characterized by three periods: 10 
a constant drying rate period, the first decreasing drying 
rate period (in which an internal wet zone with capillary 
water and an external dried zone with adsorbed water are 
formed), and a second decreasing drying period (when all 
the water is in the hygroscopic domain). Therefore, the 
following parameters were considered to calibrate the 
model. (There is also an initial warm-up period, but it 
cannot be calibrated from the final M data.) 
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The parameters B u B 2 , and B 3 are represented by the fol¬ 
lowing expression: 
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For the boundary conditions (on the external timber sur¬ 
faces), heat and vapor fluxes are assumed to be propor¬ 
tional to external heat and mass transfer coefficients, and 
the gas pressure is assumed to be equal to the kiln pressure. 


1. During the constant drying rate period, the drying rate 
depends on external factors (i.e., external heat and mass 
transfer coefficients and dry and wet bulb temperatures): 
thus a proportional factor (,S' n ) multiplying both heat and 
mass transfer coefficients is used. The adjustment used 
for the external heat transfer coefficient is shown in the 
following equation; and because the analogy between 
heat and mass was applied, 3 the same adjustment was 
used for the external mass transfer coefficient. 

Hq(stochastic) ‘^11 (8) 

2. The first decreasing drying rate period is controlled by 
heat and mass transfer throughout a dried zone that 
appears between the external surface and the still wet 
timber interior; and because heat and vapor fluxes are 
related by way of the enthalpy of vaporization, only the 
diffusive component was adjusted. Choong 11 found a lin¬ 
ear relation between the logarithm of the bound water 
diffusion coefficient and the reciprocal of the absolute 
temperature. Based on this relation, Kayihan 3 suggested 
that the random variation of the bound water diffusion 
coefficient (D b ) can be predicted with the following cor¬ 
rection factor (S b ). 

A (stochastic) ~ A e ^ ^ (9) 

3. Lastly, during the second decreasing drying rate period, 
it reduces rapidly as M tends to the equilibrium condi¬ 
tion; and this is approximately dependent on the wood- 
water activity (a). Thermodynamics indicates that the 
natural logarithm of the water activity is directly propor¬ 
tional to the difference in free energy between absorbed 
and pure water (AG). 

1 n(«) (10) 


<i - "A - l) 

(5) 

1 - Ah - A 

(6) 


(7) 


Stochastic parameters 


Because free energy change represents changes on inter¬ 
nal reversible energies, for water absorption it can be 
assumed to be approximately equal to the heat of sorp¬ 
tion. 12 Therefore, the following correction factor (S 4 ), 
proportional to the heat of sorption, was used to adjust 
the water activity. 

%,oc ha s,ic) = el5 “ ln(a)] ( 11 ) 


Because most of the model parameters produce similar ef¬ 
fects in the final result, calibration is performed by manipu¬ 
lating a small number of them but with significant and 
dissimilar effects on the final result. It is normally accepted 


Kiln simulation 

For conventional drying, timbers are placed in horizontal 
parallel layers separated by a small distance to allow pas- 
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sage of the drying air; and as the hot and dry air passes 
between timber layers, its sensible heat decreases to evapo¬ 
rate water and its relative humidity increases owing to ab¬ 
sorption of the produced vapor. These changes occur in the 
direction of the air circulation, generating a difference be¬ 
tween measured and actual temperatures at different posi¬ 
tions inside the kiln. So this difference is not mistakenly 
attributed to the external heat and mass transfer coefficient, 
the following heat and mass transfer differential equations 
are implemented to simulate the changes in the direction of 
the air flow: 


hp g v g Cp g = -2 q + 2/Cp v (7L - T s ) 

„ d JM = -2 i 

OX 


(14) 

(15) 


where the number “2” takes into account that evaporation 
occurs in both upper and lower timber surfaces. 

The kiln width in the direction of the air circulation is 
divided into a number of sectors for which the homoge¬ 
neous moisture content (M) distribution and drying air con¬ 
ditions are assumed. On each integration step, the average 
q and j fluxes on the layers are simulated simultaneously 
using several parallel stochastic models (one for each layer). 
The results are used to actualize the temperature and 
humidity profiles of the drying air. The air temperature 
(dry-bulb temperature) and humidity (or the equivalent 
wet-bulb temperature) are adjusted throughout according 
to a predetermined drying schedule, and the airflow direc¬ 
tion is reversed several times during drying. 


Table 1 . Kiln drying schedule 


Step 

Time 

(h) 

Wet bulb 
(°C) 

Dry bulb 
(°C) 

1 

6 

49 

49 

2 

24 

51 

52 

3 

24 

53 

55 

4 

24 

55 

58 

5 

24 

57 

62 

6 

24 

59 

66 

7 

24 

61 

70 

8 

24 

63 

74 

9 

90-212 

65 

78 

10 

12 

74 

78 


fore drying using a digital balance. The experimental data 
were used later to calculate M. 

Fitting to the experimental data was performed through 
numerical optimization (quasi-Newton method). In accor¬ 
dance with the most common M distribution curves found 
in the literature, 13 ' 14 experimental data were fitted on the 
basis of three key histogram characteristics: the scale 
(equivalent to the average), the shape (equivalent to the 
standard deviation), and the threshold (which quantifies the 
fact that the M distribution has a positive lower boundary). 
Four independent variables were used on the numerical 
optimization, the stochastic parameters S H , ,S' b . and S„. and 
the intrinsic dispersion (d) to fit the unknown sources of 
random behavior. 


Results and discussion 


Materials and methods 

Experimental data from six conventional dry kiln runs of 
mixed western hemlock [Tsuga heterophylla (Raf.) Sarg] 
and amabilis fir [Abies amabilis (Dougl.) Forbes], were used 
to adjust the stochastic model. 1 For each run, 160 pieces 116 
X 116mm in cross section and 2.44m long were dried to a 
target between 14% and 20% in a 7.6 m 3 conventional kiln 
installed at the University of British Columbia’s (UBC’s) 
laboratory. Timbers were arranged in eight continuous 20- 
timber layers separated by 19-mm stickers and without 
chimney space. The drying schedule was intended to be 
conservative initially yet severe enough during the latter 
stages to get a response to the timber sorting regimens 
(Table 1). Fan speed was set to 3 m/s through the timber 
tiers, but the actual speed varied between 2.95 and 3.15 m/s 
when tested. 

For determination of the M : (green) distribution, 25 mm 
thick sections were cut from each end of each green timber. 
The volume of the sections was determined by the water 
displacement method; and after oven-drying they were re- 
weighed to determine their oven-dried density (p 0 ). Green 
timber volumes were calculated from their rectangular di¬ 
mensions, and timber weights were measured after and be- 


In general, the higher the number of stochastic parameters, 
the better is the description of the final dispersion; however, 
the number of simulations increases rapidly with the param¬ 
eters, and calibration becomes ambiguous if more than one 
have a similar effect on the model’s results. Therefore, the 
optimum number of parameters is in fact the minimum 
required to provide consistent results under different ex¬ 
perimental conditions. After experimenting with several 
possibilities, it was concluded that for the conditions of this 
work (same type of wood, similar schedules, and targets 
between 14% and 20%), the final moisture content distribu¬ 
tion can be reasonably well predicted with the intrinsic 
dispersion parameter proposed in Eq. 1, so all data disper¬ 
sion introduced by several random sources could be fitted 
with only one experimental parameter independent of the 
target. Because 6 does not require more parallel simula¬ 
tions, the number of simulations is determined only by the 
Mj distribution points and the number of kiln subdivisions 
in the airflow direction. In this study three subdivisions 
(represented by three x points) were used for the kiln length 
and six subdivisions for the ;V/, distributions (represented by 
seven /W, points); thus, only 21 parallel simulations are re¬ 
quired to account for the 3x7 possible combinations. 

Two criteria were used to test the model. The first is a 
graphical comparison between experimental and simulated 
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Frequency 



M r (%) 


Fig. 1. Simulated and experimental final moisture content ( M t ) distri¬ 
bution of run 1 


Frequency 



Fig. 3. Simulated and experimental M t distribution of run 3 


Frequency Frequency 

0.35 

0.30 
0.25 
0.20 
0.15 
0.10 
0.05 
0.00 

6 9 12 15 18 21 24 27 30 33 36 

Mr (%) 

Fig. 2. Simulated and experimental M f distribution of run 2 Fig. 4. Simulated and experimental M t distribution of run 4 




Mt (%) 


histograms, which although qualitative was found to be 
more indicative of the quality of the results than other 
statistical indexes such as the mean ( X ) and standard devia¬ 
tion (o). The latter indexes were in fact used to fit the data, 
so fitted X and 6 are in all cases close to the experimental 
ones. Graphical comparisons between experimental and 
simulated M i histograms of runs 1,2,3,4,5, and 6 are shown 
in Figs. 1, 2, 3, 4, 5, and 6, respectively. It can be observed 
that an acceptable agreement exists. Exact agreement is not 
expected because the proposed model provides only the 
most probable output (the one that will be obtained after 
infinite repetitions of the same run), and the experimental 
data represent only a particular case. 

The second criterion considered is the consistency of the 
results. Even when some experimental discrepancies are 
expected due to the limited number of samples used on 
each run, the fitted coefficients for the six experimental runs 
must be similar to ensure that the model can be applied to 
predict a process with a variable target. The numerical re¬ 
sults for the fitted S H , S b , S,„ and d coefficients obtained for 


Frequency 



8 11 14 17 20 23 26 29 32 35 38 

M f (%) 

Fig. 5. Simulated and experimental M f distribution of run 5 
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Frequency 



M f (%) 


Fig. 6. Simulated and experimental M, distribution of run 6 


Table 2. Basic data and fitted stochastic parameters for the six experi¬ 
mental runs 


Run 

Time 

days 

M t 

(%) 

6 

(%) 




<5 

1 

16.6 

14.0 

4.09 

2.98 

1.47 

1.56 

0.064 

2 

15.1 

16.1 

4.65 

2.99 

1.37 

1.55 

0.062 

3 

14.1 

16.5 

5.36 

3.16 

1.49 

1.46 

0.076 

4 

12.7 

17.0 

4.63 

3.00 

1.49 

1.51 

0.063 

5 

12.8 

18.2 

5.51 

3.03 

1.34 

1.55 

0.090 

6 

11.5 

19.9 

7.00 

2.93 

1.40 

1.45 

0.076 

Mean 




3.01 

1.43 

1.51 

0.072 


M f , final moisture content; 6 , standard deviation; d, intrinsic dispersion; 
5 H , proportional factor; S b and S„, stochastic parameters 


the six runs, together with the experimental values of drying 
time, average M f , and 6 are listed in Table 2. It can be seen 
that the average S H , S b , S a , and d values are 3.01, 1.43,1.51, 
and 0.0717, respectively; and the differences among indi¬ 
vidual runs are small. For example, the 6 among each six 
parameters are 2.57%, 4.52%, and 3.08% for S H , S b , and S a , 
respectively, demonstrating the consistency of the results. 
In the case of d, the 6 of the six fitted values is higher, 
15.58%, with respect to the X , but it is still acceptable when 
taking into account that d incorporates all the random fac¬ 
tors affecting each particular run. 

Finally, the model is tested to predict target moisture 
content and its 6 using the six-run average stochastic param¬ 
eters and the six-run average M i distribution. Two com¬ 
parisons were performed between the model and the 
experimental data, namely, the average M f as a function of 
drying time (where time is only changed in step 9 of Table 
1), and 6 as a function of the average M f . The comparisons 
are shown in Figs. 7 and 8, respectively, where good agree¬ 
ment between model prediction and experimental data is 
observed. Therefore, in the range of conditions for which 
the model was adjusted, the average stochastic model 
can be used to predict conventional drying with a variable 
target. 


M t (%) 



SD (%) 



Fig. 8. Simulated and experimental standard deviation (SD) as a func¬ 
tion of M f 


Conclusions 

In this work, a stochastic model developed in previous 
works was simplified and applied to conventional drying of 
mixed western hemlock and amabilis fir timbers. Contrary 
to other, more realistic stochastic models based on a piece- 
by-piece random simulation, the proposed method takes 
samples of all the multidimensional parameter probability 
space. This new approach has two advantages: The model 
gives the same result for the same initial conditions, and the 
average and range of dispersion of the results is more or less 
well represented independent of the number of discrete 
parameter points. The stochastic model also requires imple¬ 
mentation of a deterministic timber drying model, which 
was simulated using heat and mass differential balances for 
both the timbers and the drying air. These balances are 
complicated to solve, but they incorporate flexibility to pre¬ 
dict new hypothetical situations and give information about 
the timbers internal M profiles. 
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It was found that when applied to conventional drying 
the proposed stochastic model can be considerably simpli¬ 
fied by introducing a new empirical factor that takes into 
account all sources of unknown random behavior. This new 
factor does not increase the number of simulations, so the 
model required only 21 parallel simulations to account for 
the six discrete M t subintervals and the three-length subdi¬ 
visions in the direction of the airflow. The model was tested 
with six experimental conventional runs to targets between 
14% and 20%. The results showed that the fitted param¬ 
eters and intrinsic dispersion were consistent for all runs, 
and experimental averages and standard deviations can be 
well reproduced using the six-run average parameters. 


Appendix: nomenclature 

D diffusion coefficient (m 2 /s) 

Cp heat capacity at constant pressure (J/m 3 .K) 

G free energy (J/kg) 

Hj external mass transfer coefficient (s/m) 

H q external heat transfer coefficient (W/m 2 .K) 
h separation between parallel wood layers (m) 
j external vapor flux (kg/m 2 .s) 

K intrinsic permeability of wood (m 2 ) 

K r relative permeability 

k thermal conductivity (W/m.K) 

M water moisture content (dimensionless in the model 
equations) (%) 

P pressure (Pa) 

F : ' A vapor equilibrium pressure (Pa) 

<7 external heat flux (W/m 2 ) 

R universal gas constant (J/mol.K) 

5 stochastic correction factor 

T temperature (K) 

X population mean value 

t time (s) 

v velocity (m/s) 

W molecular weight (kg/mol) 

x longitudinal in the direction of the airflow (m) 

a water activity 

d standard deviation for the unknown random 
dispersion 

p density (kg/m 2 ) 

6 moisture content standard deviation 
AH water latent heat of evaporation (J/kg) 


Subscripts 

a = relative to the air 

b = relative to the bound water diffusion coefficient 

c = relative to the capillary pressure 

eff = relative to the effective air-vapor diffusion coefficient 


f = relative to the final of drying 
g = relative to the gas (air-vapor mixture) 

H = relative to the external heat and mass transfer 
coefficient 

i = relative to the starting of drying 
1 = relative to the free water 

s = relative to the timber external surface 
v = relative to the water vapor 
a = relative to the water activity 
d = relative to the intrinsic dispersion 
°° = relative to the external drying air 
0 = relative to the oven-dried density 
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